%% First load data arrays for analysis
clear;
fnames = dir("monod*.mat");
for k=1:length(fnames)
    load(fnames(k).name);
end

clc;
close all;

%% Establish concentrations to evaluate
ivec = [1:3]'; % Only use first 4 concentrations

%% Define molar mass
molmass_succinate = 162; % Molar mass of succinate
tempmat = monod_mat_succinate_37C;
concentration = tempmat(:,1)/100*1000/molmass_succinate*1000; % Convert to mM


%% 30C succinate
tempmat = monod_mat_succinate_30C;
smat = size(tempmat);

conc_30C = tempmat(ivec,1)/100*1000/molmass_succinate*1000; % Convert to mM
grmax_30C = nanmean(tempmat(ivec,2:4),2);
grmax_err_30C = nanstd(tempmat(ivec,2:4),0,2);

%% 37C succinate
tempmat = monod_mat_succinate_37C;
smat = size(tempmat);

conc_37C = tempmat(ivec,1)/100*1000/molmass_succinate*1000; % Convert to mM
grmax_37C = nanmean(tempmat(ivec,2:4),2);
grmax_err_37C = nanstd(tempmat(ivec,2:4),0,2);

%% Now do Arrhenius analysis

conc_tot = conc_37C; % All should be the same concentration
T = [30 37]';

grmax_tot_succinate = [grmax_30C grmax_37C];
grmax_err_tot_succinate = [grmax_err_30C grmax_err_37C];

Ea_succinate = [];
Ea_err_succinate = [];

figure;
for k=1:length(ivec)
    g = grmax_tot_succinate(k,:)';
    ge = grmax_err_tot_succinate(k,:)';
    [f1,f2] = fit(1./(T+273.15), log(g),'poly1', 'Weights', (ge./g).^2);
    errorbar(1./(T+273.15), log(g), ge./g, 'ko', 'linewidth', 1);
    hold on;
    plot(1./(T+273.15),f1(1./(T+273.15)), 'k', 'linewidth', 1);
    disp('Activation Energy (KT):');
    disp(-f1.p1/298);
    
    disp('With error:');
    
    alpha = 0.90;
%     ci = confint(f1, alpha);
%     fiterr = abs(f1.p1 - ci(1,1))/298;
%     disp(fiterr);
%     
    Ea_succinate = [Ea_succinate; -f1.p1/298];
%     Ea_err_succinate = [Ea_err_succinate; fiterr];
    
    
end

set(gca, 'fontsize', 20);
set(gcf, 'Position', [0 0 400 300]);
xlabel('1/T (1/K)');
ylabel('Log (growth rate)');
xlim([1/312 1/296]);

%% Plot Ea vs. concentration (in units of kcal/mol)
figure;
errorbar(concentration, Ea_succinate/1.688, Ea_err_succinate/1.688, 'bo');
ylim([0 40]);
set(gca, 'fontsize', 20);
set(gcf, 'Position', [0 0 400 300]);
xlabel('Concentration (mM)');
ylabel('Activation energy (kcal/mol)');
box off
xlim([0 max(concentration)*1.1])

%% Compute average activation energy at saturation, with standard error
sat_index = [1:2]; % Indexing for saturating concentrations

Ea_mean = mean(Ea_succinate(sat_index))/1.688
Ea_mean_err = sqrt(sum(Ea_err_succinate(sat_index).^2)/length(sat_index))/1.688
% Ea_se = Ea_mean_err/sqrt(length(sat_index));









